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Additional Keplerian Signals in the HARPS data for 
Gliese 667C from a Bayesian Re-analysis 
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ABSTRACT 

A re-analysis of Gliese 667C HARPS precision radial velocity data was carried out with 
a Bayesian multi-planet Kepler periodogram (from to 7 planets) based on a fusion 
Markov chain Monte Carlo algorithm. The most probable number of signals detected 
is 6 with a Bayesian false alarm probability of 0.012. The residuals are shown to be 
consistent with white noise. The 6 signals detected include two previously reported 
with periods of 7.198 (b) and 28.14 (c) days, plus additional periods of 30.82 (d), 38.82 
(e), 53.22, and 91.3 (f) days. The 53 day signal is probably the second harmonic of 
the stellar rotation period and is likely the result of surface activity. The existence of 
the additonal Keplerian signals suggest the possibilty of further planets, two of which 
(d and e) could join Gl 667Gc in the central region of the habitable zone. N-body 
simulations are required to determine which of these signals are consistent with a 
stable planetary system. Msinz values corresponding to signals b, c, d, e, and f are ~ 
5.4, 4.8, 3.1, 2.4, and 5.4 M^, respectively. 

Key words: stars: planetary systems; methods: statistical; methods: data analysis; 
techniques: radial velocities. 



1 INTRODUCTION 

The HARPS spectrograph mounted on the ESO/3.6-m La 
Silla telescope is yielding a rich bount y of information re - 
garding planets around M dwarfs fee.. iBonfils et ailbOllT ). 
There is a lot of interest in searches for planets about low 
mass M dwarfs. Firstly, a planet of given mass and orbital 
separation induces a larger stellar radial velocity variation 
around a lower mass star. Secondly, the low luminosity of 
M dwarfs moves their habitable zone much nearer to the 
star. For these reasons a habitable zone (HZ) planet around 
a 0.3 Mq M dwarf produces a 7 times larger radial velocity 
wobble than the sam e planet orbiting a solar-mass G dwarf 
(jPelfosse et al.ll2012l ). 

One M dwarf of particular interest is Gliese 667C (Gl 
667C, also known as GJ 667C), an isolated component of a 
hierarchical triple system. The two others components, Gl 
667AB, are a closer couple of K dwarfs. Gl 667AB has a 
semi-major axis of 1.82 AU (period of 42.15 years) and a to- 
tal m ass of 1.27 Mq (dynamically determined, Sodcrhiclm 
1 19991 ). Gl 667C is the lightest component with a mass of 
0.310 ± 0.019 M0 (|Anglada-Escude et al.ll2012al '). R is at a 
projected distance of 32.4" from Gl 667AB, giving an ex- 
pected semi-major axis of ~ 300 AU (for a distance of 7.23 
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pc and a fac tor of 1.26 between expec ted and projected semi- 
major axis, iFischer and Marcvlfigg^ l. 



In 2011. lBonfils et al.l |2oT3) reported a super-earth Gl 
667Cb (Msini = 5.9 M0) with a period of 7.2 d and 
evidence for two other interesting periods of 90 and 28 
d. The possibility of a 28 d period planet was particu- 
larly interest ing because it would fall in th e HZ. Two re- 
cent papers (|Anglada-Escude et al.l l2012al & iDelfosse et al.l 
[2013' ), have confirmed planet Gl 667Cb and a 28 d pe- 
riod planet Gl 667Gc (Msini = 4.3 M0) in the HZ. The 
lAnglada-Escude et al.l l|2012ar ) results are not full y inde - 
pendent as they depend on the 143 iBonfils et all (|201lD 
RV observations which they reduced using the HARPS- 
TERRA (Template En hanced Radial velocity Re-an alysis 
Application) algorithm (|Anglada-Escude et al.ll2012bl ). sup- 
ported by obs e rvatio ns from two other spectrographs. The 
iDelfosse et"aLll|2012h data includes an additional 36 HARPS 
radial velocity measurements. 

The excitement generated by this and many other exo- 
planetary discoveries has spurred a significant effort to im- 
prove th^_stati^tical__to^ 



(e.g.. lLoredo fc Ghernofj|2003l.lLoredoll2004ICummindl2004 



Gregory 2005a fc b. Ford 2005 fc 2006. iFord fc Gregory! 
2007I. iTuomi fc Kotirai^ '20091. iDawson fc Fabrvckvl 120101 



Gumming fc Dragomir 2010. ) . Much of this work has high- 
lighted a Bayesian MCMC approach as a way to better un- 
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derstand parameter uncertainties and degeneracies and to 
compute model probabilities. 

Gregory (2011, 2011b, 2012) developed a Bayesian fu- 
sion MCMC algorithm that incorporates parallel tempering 
(PT), simulated annealing and a genetic crossover operation 
to facilitate the detection of a global minimum in x^- This 
enables the efficient exploration of a large model parameter 
space starting from a random location. When implemented 
with a multi-planet Kepler model Q, it is able to identify 
any significant periodic signal component in the data that 
satisfies Kepler's laws and is able to function as a multi- 
planet Kepler periodogram 0. The fusion MCMC algorithm 
includes an innovative adaptive control system that auto- 
mates the selection of efficient parameter proposal distribu- 
tions even if the parameters arc highly corre lated . One appli- 
cation of the algorithm ([Greeorv fc Fischeili2010|) confirmed 
the existence of a disputed second planet I Fischer et al.l 
I2OO2I ') in 47 Ursae Majoris (47 UMa) and provided orbital 
constraints on a possible additional long period planet with 
a period ~ lOOOOd. 

This pape r presents a Bayes ian re-analysis of the 
HARPS data (jPelfosse et all l2012h for the star Gl 667C 
using our fusion MCMC based multi-planet Kepler peri- 
odogram. Section[2]provides a brief overview of our Bayesian 
approa ch and outlines the adaptive fusi on MCMC algo- 
rithm (ICregorvllibllbl and [Cregorvl 120121 '). Section [3] gives 
the model equations and priors. Sections [4] and [5] present 
the parameter estimation and model selection results. The 
final two sections are devoted to discussion and conclusions. 



2 THE ADAPTIVE FUSION MCMC 

The adaptive fusion MCMC (FMCMC) is a very general 
Bayesian nonlinear model fitting program. After specify- 
ing the model, Mi, the data, D, and priors, /, Bayes theo- 
rem dictates the target joint probability distribution for the 
model parameters which is given by 



p{X\D,M„I)^C p{X\M„I)xp{D\M„X,I). 



(1) 



where C is the normalization constant and X represent the 
set of model parameters. The first term on the RHS of the 
equation, p{X\Mi, I), is the prior probability distribution of 
X, prior to the consideration of the current data D. The 
second term, p(D\X , Mi, I), is called the likelihood and it is 
the probability that we would have obtained the measured 
data D for this particular choice of parameter vector X, 
model Mi, and prior information I. At the very least, the 
prior information, /, must specify the class of alternative 
models (hypotheses) being considered (hypothesis space of 
interest) and the relationship between the models and the 
data (how to compute the likelihood). In some simple cases 
the log of the likelihood is simply proportional to the familiar 



^ For multiple planet models, there is no analytic expression for 
the exact radial velocity perturbation. In many cases, the radial 
velocity perturbation can be well modeled as the sum of multiple 
independent Keplerian orbits which is what has been assumed in 
this paper. 

^ Following on from the pion eering work on Bayesian peri- 
odograms bv ljavnej l ll987t) and iBretthorsj lll988l ') 



statistic. For further details of the likelihood function for 
this type of problem see Gregory C2Q0 5b). 

To compute the marginals for any subset of the 
parameters it is necessary to integrate the joint probability 
distribution over the remaining parameters. For example, 
the marginal probability distribution (a probability density 
function) of the orbital period in a one planet radial velocity 
model fit is given by 



p{P\D,M^,I) = 



dK dV de dx 



X I dcj J dslope J ds 
X p{P, K, V, e, X, t^, slope, s\D, Mi, I) 
(X p(P|Mi,/) dK - - I ds 



X p{K, V, e, X, ^, slope, s\Mi, I) 

X p{D\Mi,P,K,V,e,x,u},slope,s,I), 



(2) 



where P, K,V, e,x,'jJ are the radial velocity model parame- 
ters and s is an extra noise parameter which is discussed in 
Section [3l slope is a parameter that accounts for the long 
period orbital motion induced in Gl 667C by Gl 667AB. 
p{P\Mi, I) is the prior for the orbital period parameter, 
p{K,V, e,x,^, slope, s\Mi, I) is the joint prior for the other 
parameters, and p{D\Mi, P, K,V, e,x,^-^, slope, s, I) is the 
likelihood. For a seven planet model fit we need to inte- 
grate over 37 parameters to obtain the marginal probability 
distribution for each of the seven period parameters. Inte- 
gration is more difficult than maximization, however, the 
Bayesian solution provides the most accurate information 
about the parameter errors and correlations without the 
need for any additional calculations, i.e., Monte Carlo sim- 
ulations. Bayesian model selection requires integrating over 
all the model parameters. 

In high dimensions, the principle tool for carrying out 
the integrals is Markov chain Monte Carlo based on the 
Metropolis algorithm. The greater efficiency of an MCMC 
stems from its ability, after an initial burn-in period, to gen- 
erate samples in parameter space in direct proportion to the 
joint target probability distribution. In contrast, straight 
Monte Carlo integration randomly samples the parameter 
space and wastes most of its time sampling regions of very 
low probability. 

MCMC algorithms avoid the requirement for com- 
pletely independent samples, by constructing a kind of ran- 
dom walk in the model parameter space such that the num- 
ber of samples in a particular region of this space is pro- 
portional to a target posterior density for that region. The 
random walk is accomplished using a Markov chain, whereby 
the new sample, Xt+i, depends on previous sample Xt ac- 
cording to a time independent entity called the transition 
kernel, p{Xt+i\Xt)- The remarkable property of p{Xt+i\Xt) 
is that after an initial burn-in period (which is discarded) 
it generates samples of X with a probability density pro- 
portional to t he desired poste rior p{X\D, Mi, I) (e.g., see 
Chapter 12 of lCregorvl (|2005al ) for details). 
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The joint posterior probability distribution in model pa- 
rameter space is typically highly multi-modal for exoplanet 
radial velocity (RV) analysis. The Metropolis algorithm can 
become stuck in the vicinity of a local probability maxi- 
mum. To avoid this fusion MCMC (FMCMC) incorporates 
three other algorithms each of which is designed to facilitate 
the detection of a global minimum in (or a maximum in 
probability). They are parallel tempering (PT) (lGeverlll99ll 
and re-invented bv lHukushima fc Nemotoll 19961 ). simulated 
annealing and a genetic crossover operation. All three are 
implement in each FMCMC run. The combination greatly 
facilitates the identification of a global maximum in prob- 
ability. This was made possible through the development 
of an adaptive control system which has been described in 
detail most recently by Gregory (2011b, 2012). Further re- 
finements of the control system are ongoing. 

At each iteration of the FMCMC, a single joint pro- 
posal is made to jump to a new location in parameter space 
from the current location. The key to an efficient MCMC 
is an efficient method of proposing new jumps especially 
when there are correlations between the parameters. This 
is further complicated in PT because multiple MCMC are 
executed in parallel each exploring a different probability 
distribution. In FMCMC this difficult step is automated by 
the control system saving a great deal of time and effort. 



3 MODELS AND PRIORS 

In this section we describe the model fitting equations and 
the selection of priors for the model parameters. We have 
investigated the Gl 667C data using models ranging from 
to 7 planets. For a one planet model the predicted radial 
velocity is given by 

v{U) = 1/ + K[cos{e{U + xP) + ^} + ecoso;], (3) 

and involves the 6 unknown parameters 

V — a constant velocity. 

K — velocity semi-amplitude. 

P — the orbital period. 

e — the orbital eccentricity. 

oj = the longitude of periastron. 

X = the fraction of an orbit, prior to the start of data 
taking, that periastron occurred at. Thus, xP = the number 
of days prior to ti = that the star was at periastron, for 
an orbital period of P days. 

Qiti 4- xP^ ~ the true anomaly, the angle of the star in 
its orbit relative to periastron at time ti. 

We utilize this form of the equation because we obtain 
the dependence of Q on ti by solving the conservation of 
angular momentum equation 

dB_ _ 27r[l + ecosg(t.+xP)l' ^ , . 

dt P(l -62)3/2 

Our algorithm is implemented in Mathematica and it proves 
faster for Mathematica to solve this differential equation 
than solve the equations relating the true anomaly to the 
mean anomaly via the eccentric anomaly. Mathematica gen- 
erates an accurate interpolating function between t and 9 
so the differential equation does not need to be solved sepa- 
rately for each ti. Evaluating the interpolating function for 



each ti is very fast compared to solving the differential equa- 
tion. Details on how equation |4] is implemented and the ac- 
curacy of the interpolation as a func t ion of eccentricity are 
given in the Appendix A of iGregor vl ll2011bl) 

As described in more detail in iGregorvllioOTl . we em- 
ployed a re-parameterization of x s-^d uj to improv e the 
MCMC convergence speed motivated by the work of iFordI 
l|2006l ). The two new parameters are ip = 2-rx + ^ S'ld 
(j> = '2-KX — ^- Parameter is well determined for all ec- 
centricities. Although (j) is not well determined for low ec- 
centricities, it is at least orthogonal to the 'i/' parameter. 
We use a uniform prior for tp in the interval to 4t7 and 
uniform prior for in the interval ~2n to +2n. This in- 
sures that a prior that is wraparound continuous in (x, w) 
maps into a wraparound continuous distribution in {il),(j>). 
To account for the Jacobian of this re-parameterization it 
is necessary to multiply the Bayesian integrals by a factor 
of (47r) where nplan = the number of planets in the 
model. In this work we utilized the orthogonal combination 
(i/', (p) as we ll as the FMCMC correlated proposal scheme 
described in iGregorvl (|2012l ). 

In a Bayesian analysis we need to specify a suitable prior 
for each parameter. These are tabulated in Table [T] For the 
current problem, the prior given in equation[2]is the product 
of the individual parameter priors. De tailed argumen ts for 
the choice of each prior were given in IGregorvl l|2007h and 
iGregorv fc FischeJ (|2010l ). 

All of the models considered in this paper incorporate 
an extra noise parameter, s, that can allow for any addi- 
tional noise beyond the known measurement uncertainties. 
We adopt an independent Gaussian distribution with a vari- 
ance where is becomes another parameter in the model 
fit. Thus, the combination of the known errors and extra 
noise has a Gaussian distribution with variance = + s^, 
where ai is the known measurement uncertainty for i"' data 
point. In general, nature is more complicated than our model 
and known measurement errors. Marginalizing s has the de- 
sirable effect of treating anything in the data that can't be 
explained by the model and known measurement errors as 
noise, leading to conservative estimat es of orbi t al para me- 
tersB See Sections 9.2.3 and 9.2.4 of lGregorvl (jlooHi) for 
a tutorial demonstration of this point. If there is no extra 
noise then the posterior probability distribution for s will 
peak at s = 0. The upper limit on s was set equal to A'max. 
We employed a modified scale invariant prior for s with a 
knee, so = 1 m s~^. 



3 There is an additional benefit for incorporating the extra noise 
parameter s. When the of the fit is very large, the Bayesian 
Markov chain automatically inflates s to include anything in the 
data that cannot be accounted for by the model with the current 
set of parameters and the known measurement errors. This results 
in a smoothing out of the detailed structure in the surface 
and allows the Markov chain to explore the large scale structure 
in parameter space more quickly. This is similar to simulated 
annealing, but does not require choosing a cooling scheme. The 
chain begins to decrease the value of the extra noise as it settles 
in near the best -fit parameters. This i s demonstrated in Fig. 2 of 
IGregorvl l l2011bl ') and in Section 2.1 of lGregorvl | |2012|) . 
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Table 1. Prior parameter probability distributions. 



Parameter 



Prior 



Lower bound Upper bound 



Orbital frequency 

Velocity Ki 
(m s-i) 

V (m s-l) 

Si Eccentricity 

X orbit fraction 



LJi Longitude of 
periastron 



slope 

m/s/yr 

s Extra noise (m s~^) 



p(ln/i,ln/2,---ln/„|M„,J) = 
(n =number of planets) 



Modified scale invariant 



(y+Ko)-i 



1/3 



Uniform 

Ecc. noise bias correction filter 

Uniform 

Uniform 

Uniform^ 



1/0.5 d 1/95 yr_^ 



^ 1/3 



(Ko = 1) /fmax ) 



A-max = 2129 



0.99 
1 

27r 



(so = 1) Kn 



Gl 667C is part of a triple star system with Gl 667AB, a much closer binary. We adopted at an upper 
limit of 95 yr by sotting the gravitational pull on the planet due to Gl 667AB = 1% of the pull from Gl 
667C. 

^ Since the prior lower limits for K and s include zero, we used a modified scale invariant prior of the form 

p(X\M,I) = 2— (5) 

X + Xo ln(l + 2^) 

For X ■C Xo, p{X\M, I) behaves like a uniform prior and for X S> it behaves like a scale invariant prior. 
The In (l + ^^^J^) term in the denominator ensures that the prior is normalized in the interval to Xmax- 
Since this parameter is common to all models including the planet model, the exact upper and lower 
bounds are not critical for model selection. The range simply needs to be large enough so as not to effect 
the parameter estimates. 



3.1 Eccentricity bias 

In iGregorv fc FischeJ (|2010t ). the velocities of model fit 
residuals were randomized in multiple trials and processed 
using a one planet version of the FMCMC Kepler peri- 
odogram. In this situation periodogram probability peaks 
are purely the result of the effective noise. The orbits corre- 
sponding to these noise induced periodogram peaks exhib- 
ited a well defined statistical bias Q towards high eccentric- 
ity. We offered the following explanation for this effect. To 
mimic a circular velocity orbit the noise points need to be 
correlated over a larger fraction of the orbit than they do 
to mimic a highly eccentric orbit. For this reason it is more 
likely that noise will give rise to spurious highly eccentric 
orbits than low eccentricity orbits. 



* The bias found using multiple sets of randomized residuals from 
a 5 planet fit to 55 Cancri combined Lick and Keck data agreed 
closely with the bias found for multiple sets of randomized resid- 
uals from both 2 and 3 planet fits to 47 UMa Lick data. 



iGreeorv fc FischeJ (|201(]| ) characterized this eccentric- 
ity bias and designed a correction filter that can be used as 
an alternate prior for eccentricity to enhance the detection 
of planetary orbits of low or moderate eccentricity. On the 
basis of our understanding of the mechanism underlying the 
eccentricity bias, we expect the eccentricity prior filter to 
be generally applicable to searches for low amplitude orbital 
signals in precision radial velocity data sets. The probability 
density function for this filter is given by 

pdf{e) = 1.3889 - 1.5212e^ + 0.53944e^ 

-1.6605(e- 0.24821)*. (6) 

Fig. 11 of iGregorvl (|2011bl ) demonstrates that the effect 
of this noise eccentricity bias correction filter on the final 
marginal eccentricity distributions for the low and moder- 
ate eccentricity orbits of Gl 581b, c, fc d is insign ificant. 

In a related study, IShen and Turneil (|2008h explored 
least-x^ Keplerian fits to synthetic radial velocity data sets. 
They found that the best fit eccentricities for low signal-to- 
noise ratio K/g ^ 3 and moderate number of observations 
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-1000 -500 500 1000 

Julian day number (-2,454,504.8055) 

Figure 1. The HARPS data for Gl 667C after correcting for the 
secular acceleration. 



Nobs ^ 60, were systematically biased to higher values, lead- 
ing to a suppr ession of the number of nearly circular orbits. 
More recently, IZakamska et al.l (|201ll ') found that eccentric- 
ities of planets on nearly circular orbits are preferentially 
overestimated, with typical bieis of one to two times the me- 
dian eccentricity uncer tainty in a survey, e.g., 0.04 in the 
Butler et al. catalogue (|Butler et al. I I2OO6I '). When perform- 
ing population analysis, they recommend using the mode of 
the marginalized posterior eccentricity distribution to min- 
imize potential biases. 

In the analysis of the Gl 667C data we used the eccen- 
tricity noise bias correction filter as the eccentricity prior on 
fits of all the models. 



4 ANALYSIS OF THE HARPS DATA 

For all the analysis we used the HARPS data given by 
iDelfosse et al] (|2012l ). We subtracted a mean velocity of 
6.5 474477 km s~^ and co nverted to units of m s~^. Follow- 
ing iDelfosse et al] l|2012l ) we also subtracted a component 
due to the secular acceleration of 0.21 m s~^ yr"^ arising 
from the proper motion. This small correction resulted in a 
change in radial velocity over the span of the data amount- 
ing to 1.527 m s~^. Fig.[l]shows the corrected HARPS data 
for Gl 667C which exhibits a large positive slope indicative 
of a period much longer than the data. A likely explana- 
tion is the orbital motion of Gl 667C relative to the center 
of mass of the Gl 667ABC system. The Gl 667ABC sys- 
tem parameters suggested an orbital period of ~ 3900 yr. 
For a circular orbit this implied a maximum expected K 
value for Gl 667C of ~ 1850 m s~^ and a maximum rate 
of change in radial velocity of ~ 3 m s~^ yr~^. Multiplying 
by the data duration of 7.27 yr yields a maximum expected 
velocity change of 20 m s~^, which is comparable to what 
is observed. Our zero reference time is the mean time of 
the HARPS observations which corre sponds to a Jul i an da y 
number = 2, 454, 504.8055. Following lOelfosse et al.1 l|2012ll . 
we deal with the long period velocity change by including a 
linear drift which we refer to as a slope parameter in all our 
models. 

The models considered range from a zero planet model 
to a seven planet model. As mentioned all models include a 
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Figure 2. The upper panel is a plot of the Logio [Prior X Like- 
lihood] versus iteration for a 2 planet Kepler periodogram of the 
HARPS data. The lower shows the values of the two unknown 
period parameters versus iteration number. The two starting pe- 
riods of 7.2 and 15 d are shown on the left hand side of the plot 
at a negative iteration number. 



slope parameter and extra noise (s) parameter. The current 
section deals with model parameter estimation while Sec- 
tion [5] deals with model selection. The next 6 subsections 
show the FMCMC based Kepler periodograms for the 2-7 
planet cases. We don't show Kepler periodogram for the one 
planet model as its parameters are well understood. 



4.1 Two planet model 

Five 2 planet Kepler periodograms were computed. In all 
cases signals were detected at 7.2 and 28 d but the 28 d 
signal was never the dominant second peak. Other periods 
which occured in individual periodograms included 106, 184, 
249, 383 d. Fig. [2] shows the results for the periodogram 
which achieved the highest peak Logio [Prior x Likelihood] . 
The upper plot shows the Logio [Prior x Likelihood] versus 
iteration. The lower plot shows the two period parameters 
versus iteration which shows a stable 7.2 d signal while the 
second period transitions over many peaks thanks to the 
parallel tempering feature of the FMCMC algorithm. The 
relative peak probabilities of the second period options is 
illustrated in Fig. [3] which shows a plot of the two period 
parameter values versus a normalized value of Logio [Prior x 
Likelihood]. Fig. U shows a plot of the eccentricity parameter 
values versus period. Only the 7.2 and 28.1 d periods are 
consistent with low eccentricity orbits. The median value of 
the extra noise parameter s = 1.8 m s~^. 



4.2 Three planet model 

Six 3 planet Kepler periodograms were computed for the 
data. In all cases signals were detected at 7.2 and 28.1 d. 
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28 106 

Period (d) 



186 249 374 



Figure 3. A plot of the 2 period FMCMC parameter samples 
versus a normalized value of Logio [Prior X Likelihood] for the 2 
planet Kepler periodogram. 




Figure 6. A plot of eccentricity versus period for the FMCMC 
parameter samples from the 4 planet Kepler periodogram.. 



28 106 186 249 374 7451000 

Periods 




Figure 4. A plot of eccentricity versus period for the FMCMC 
parameter samples from the 2 planet Kepler periodogram. 



Figure 7. A plot of the five period parameter values versus a 
normalized value of Logio [Prior X Likelihood] for the 5 planet 
Kepler periodogram. 



Other periods which occured in individual periodograms 
were 38.8, 106, 184, 368 d. In three cases the third period 
was 184 d. Fig. [5] shows a plot of eccentricity versus pe- 
riod for the periodogram which achieved the highest peak 
Logio [Prior x Likelihood]. The 7.2 and 28.1 d signals are 
consistent with low eccentricity orbits. The third period in- 
variably had a high eccentricity as is the case for the 184 
d period in the figure. The median value of the extra noise 
parameter s = 1.4 m s~^. 



4.3 Four planet model 

Three 4 planet Kepler periodograms were computed for the 
data. In all cases signals were detected at 7.2, 28.1, 106, & 
190 days. Additional periods that were detected in different 
runs were 38.8 and 91d. Fig. |6] shows a plot of eccentricity 
versus period for the periodogram which achieved the high- 
est peak Logio [Prior x Likelihood] . The median value of the 
extra noise parameter s = 1.0 m s~^. 



0.8 ■ 

^ 0.6 ■ 
o 

0) 

o 0.4 ■ 
LU 

0.2 ■ 

0.0 1 , , , , , . : 

1 3 7.2 28 100 184 300 1000 

Periods 

Figure 5. A plot of eccentricity versus period for the FMCMC 
parameter samples from the 3 planet Kepler periodogram. 



4.4 Five planet model 

Four 5 planet Kepler periodograms were computed for the 
data. In all cases signals were detected at 7.2, 28.1, 53.2, & 91 
d. Other periods which occured in individual periodograms 
included 30.8, 38.8, 87, 106, 128, 184 d. In some cases the 
fifth period exhibited a variety of periods. Fig. [7] shows the 
periodogram which achieved the highest peak Logio [Prior x 
Likelihood] . In this case the extra period exhibits a dominant 
peak at 86.6 d and a weaker peak at 30.8 d. Fig. |8]is a plot of 
eccentricity versus period. The 86.6 d signal exhibits a high 
eccentricity while the 30.8 d signal has a low eccentricity. 
The median value of the extra noise parameter s = 0.79 m 
s-i. 
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Figure 8. A plot of eccentricity versus period for the FMCMC 
parameter samples from the 5 planet Kepler periodogram.. 



4.5 Six planet model 

We next computed a 6 planet Kepler periodogram analysis 
of the data. An initial blind search with starting periods 
of 5, 10, 15, 20, 50, 100 d yielded 5 well defined periods 
at 7.2, 28.1, 30.8, 53.2, 91.3 d. The remaining period was 
split between a high eccentricity 14.4 d period and a low 
eccentricity 35 d period. A second run starting from the 5 
well defined periods plus the 35 d low eccentricity option is 
shown in Fig. [9] After a temporary stay at the 35 d period, 
the black trace transitions to a stable 38. 8 d peak 0. As 
discussed in Section [5] on model selection, the Bayes factor 
favors the 6 planet model by a factor of 137 compared to 
the next leading contender. The median value of the extra 
noise parameter s — 0.54 m s~^. 

Fig. [9] shows a plot of the 6 period parameter values (in- 
cluding burn-in points to show the 35 d signal) versus a nor- 
malized value of Logio [Prior x Likelihood] for the 6 planet 
Kepler periodogram. The fourth period, shown in black, ex- 
hibits two peaks, a weak one at 35 d and the other at 38.8 
d. The 35 d period coincides with a one year alias of the 
stronger 38.8 d period. The spectral window function for the 
HARPS data exhibits two peaks, 1 d and one year. We can 
gain further insight into the relationship of the 35 and 38.8 
d periods from a 6 planet Kepler periodogram of a subset of 
the HARPS data, i.e., the first 143 data points. Again dur- 
ing the burn-in, the fourth period locks on to the 35 d peak 
before transitioning to the 38.8 d peak. Fig. [11] shows the 
6 period parameter values (including burn-in points) versus 
a normalized value of Logio [Prior x Likelihood] for the 143 
d sub-sample. In this case, the alias at 35 d actually has a 
higher peak probability density by a factor of ~ 10, even 
though the much larger number of samples for the 38.8 d 
peak indicates that there is more probable associated with 
the 38.8 d peak. A narrow peak in the joint parameter space 
of high probability density can contain much less total prob- 
ability than a broader region of lower probability density. 
Note that the other one year alias of the 38.8 d period at 



^ Another 6 planet Kepler periodogram was computed starting 
from the best 6 planet parameters set with the exception that 
53.2 d period which was replaced by twice this period which cor- 
respon ds roughly with the assumed rotation period of the star of 
105 d jPelfosse et al.ll2012l . based on a periodogram of a stellar 
activity diagnostic). This solution with periods of 7.2, 28.1, 30.8, 
38.8, 91 and 106.5 d had a peak probability 2600 times lower. 




0.2 0.3 

Iterations (« 10^) 




0.2 0.3 
Iterations (x 10°) 

Figure 9. The upper panel is a plot of the Logio [Prior x Like- 
lihood] versus iteration for the six planet fit of the HARPS data. 
The lower shows the values of the six unknown period parameters 
versus iteration number. The six starting periods are shown on 
the left hand side of the plot at a negative iteration number. 
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Period (d) 



300 



Figure 10. Full HARPS data set: A plot of the 6 period parame- 
ter values versus a normalized value of Logio [Prior X Likelihood] 
for the 6 planet Kepler periodogram (includes burn-in points to 
show the 35 d signal). The fourth period, shown in black, exhibits 
two peaks one at 35 d and the other at 38.8 d. The 35 d period 
is a one year alias of the stronger 38.8 d period. 



43.4 d coincides with the location of a very weak feature. 
Comparing Fig's [TO] and 1111 it is clear that the addition 
of the extra data p oints in the latest HARPS sample of 
[Pelfosse et al.l l|2012i ) has suppressed the alias at 35 d by a 
factor of ~ 10^ 

It was not until the 6 planet model that the 30.8 and 
38.8 d signals became the dominant third and fourth periods. 
Both are consistent with low eccentricity orbits as shown in 
Fig. [12] which is a plot of eccentricity versus period for the 
post-burn-in samples for the full HARPS data set. 

The closeness of the 28.1 and 30.8 day periods raises 
questions about a possible relationship between these two 
signals. Could the 30.8 d signal be an alias of the 28.1 d 
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3. 7.2 16 28 38.8 54 91 

Period (d) 

Figure 11. Subset of HARPS data (first 143 points): A plot 
of tile 6 period parameter vaiues versus a normalized value of 
Logio [Prior X Likelihood] for the 6 planet Kepler periodogram 
(includes burn- in points to show the 35 d signal). The fourth 
period, shown in black, exhibits two peaks one at 35 d and the 
other at 38.8 d. For the partial data set the 35 d alias has a higher 
peak probability density. The other one year alias at 43.4 d is just 
discernible. 




28.1 38.8 53.2 91.3 

Periods 



Figure 12. Full HARPS data set: A plot of eccentricity versus 
period for the 6 planet Kepler periodogram for the post-burn-in 
FMCMC samples. 



signal? An alias results from a convolution of the spectral 
window function of the data sample times with the spec- 
trum of the real signal. If the real signal is removed then an 
alias of that signal should not be present in the residuals. 
In our case the multi-planet model is requiring both to be 
present. Nevertheless, we transformed the marginal distri- 
bution of the second period, P2 (28.14 d signal), by the rela- 
tionship Paiias = 1/(1/^2 — 1/365.25) to obtain the pridicted 
marginal distribution the one year alias. This is compared 
it to the marginal for P3 (30.82 d signal) in the top panel 
of Fig. 1131 Clearly, the predicted alias (shown dashed) does 
not overlap the marginal for P3. For comparison, the bottom 
panel shows the predicted marginal distribution of the lower 
one year alias (dashed) of the 38.8 d signal to the marginal 
distribution obtained for the 35 d signal (solid) which was 
observed during the burn-in phase of the full HARPS data 
set. In this case the distributions coincide. 

Fig. [Til shows a plot of a subset of the FMCMC param- 
eter marginal distributions for the 6 planet FMCMC fit of 
the data. The eccentricities of the 4 lowest periods are con- 
sistent with low eccentricity orbits while the 53.2 and 91.3 
d periods show peak eccentricities of 0.41 and 0.36, respec- 




30.6 30.7 30.8 
P (d) 




34.8 34.9 35.0 35.1 35.2 35.3 35.4 35.5 

P (d) 

Figure 13. The upper panel shows a comparison of the predicted 
marginal distribution (dashed) of the closest one year alias of the 
28.14 d signal to the marginal distribution for the 30.82 d signal 
(solid). For comparison, the bottom panel shows the predicted 
marginal distribution of the lower one year alias (dashed) of the 
38.8 d signal to the marginal distribution obtained for the 35 d 
signal (solid) which was observed during the burn-in phase of the 
full HARPS data set. 



lively. The median value of extra noise parameter s = 0.54 
m s^^ . 

Phase plots for the 6 planet model are shown in Fig. 1151 
The top left panel shows the data and model fit versus 7.2 
d period phase after removing the effects of the five other 
orbital periods plus V and slope parameter. The FMCMC 
output for each iteration is a vector of the 6 planet orbital 
parameter set plus V , slope, and the extra noise parameters. 
The 7.2 d period phase plot is constructed from a sample of 
FMCMC iterations (typically 300) and for each iteration we 
compute the predicted velocity points for that realization 
of the 5 planet plus V and slope parameters. We then con- 
struct the mean of these model predictions and subtract the 
mean prediction from the data points. These residuals for 
the set of observation times are converted to residuals ver- 
sus phase using the mode of the marginal distribution for the 
7.2 d period parameter. A period phase model velocity fit 
is then computed at 100 phase points for each realization of 
the 7.2 d planet parameter set obtained from the same sam- 
ple of 300 FMCMC iterations. At each of these 100 phase 
points we construct the mean model velocity fit and mean 
±1 standard deviation. The upper and lower solid curves in 
Fig. mi are the mean FMCMC model fit ±1 standard devi- 
ation. Thus, 68.3% of the FMCMC model fits fall between 
these two curves. 

The next five panels correspond to phase plot for the 
other 5 periods. In each panel the quoted period is the mode 
of the marginal distribution. The bottom panel is for the 
slope parameter after removing the effects of the 6 orbital 
periods plus V. 

Fig. [16] shows a generalized Lomb-Scargle (GLS) peri- 
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Figure 14. A plot of a subset of the FMCMC parameter marginal 
distributions for the 6 planet Kepler periodogram. 



odoeram llZechmeister fc Kiirste j [ioogl 'l for the maximum 
a posteriori (MAP) parameter values of the 6 planet fit 
residuals. The GLS allows for a floating offset and weights. 
The dashed horizontal lines correspond to peak periodogram 
levels for which the false alarm probability (FAP) would 
= 0.1 & 0.010. There is no evidence of any significant peaks 
or red noise in these residuals. 

Another way of exploring the residuals is to compute 
the autocorrelation function, p{j). Fig. ll7l shows p[j) of the 
residuals for the 6, 3 and 1 planet fits computed from equa- 
tion (O. 



P{j) = 



E 



overlap [(^' i^^+J-^)] 



overlap 



-/overlap 



{xi+j ~ a;) 2 



(7) 



where Xi is the i^^ residual, j is the lag and x is the mean 
of the samples in the overlap region. Because the data are 
not uniformly sampled, for each lag all sample pairs that 
differed in time by this lag ±0.1 d were utilized. The bottom 

^ The interesting region of power is where the frequentist p- value 
is small (<IC 1). In this regio n the FAP is given approxim ately by 
FAPRi M * (1 - z)('V-3)/2 ||Zechmeister fc Kiirste3l200gh . where 
z = maximum periodogram power, M = the number of indepen- 
dent frequencies and A'^ is the number of data points. Clearly, 
the FAP value for the a ctual highe s t pea k in the periodogram 
is much larger than 0.1. ICumming l |2004 ) recommends setting 
M = Af/6f, where A/ is the frequency range examined /max, 
and 5f = the resolution of the periodogram Ri 1/T where T is 
the duration of the data. 




Figure 15. The top left panel shows the data and model fit versus 
7.2 d period phase after removing the effects of the 5 other orbital 
periods plus V and slope parameters. The upper and lower curves 
are the mean FMCMC model fit ±1 standard deviation. The next 
five panels correspond to phase plot for the other five periods. The 
bottom panel is for the slope parameter after removing the effects 
of the 6 orbital periods plus V. 




1.0 

frequency (1/d) 

Figure 16. A periodogram of the 6 planet fit residuals. 



right panel of Fig. \T7\ shows a plot of the number of such 
sample pairs as a function of the lag. Clearly for large lags 
the uncertainty in computed p{j) is expected to be larger. 

Clearly, the 6 planet model autocorrelation function is 
consistent with white noise. The thick solid line in the 1 
planet residuals panel is the average autocorrelation func- 
tion generated from 400 simulated data sets of a 5 planet 
model (28.1, 30.8, 38.8, 53.2, and 91.4 d periods) together 
with the measurement errors. The 5 planet model param- 
eter values are the MAP values derived from the 6 planet 
Kepler periodogram of the real data. The good agreement 
between the simulation and actual p{j) argues that in the 
case of Gl 667C any colored noise evident in the 1 planet fit 
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Figure 17. The top panel shows the autocorrelation function, 
p{j), of the HARPS 6 planet fit residuals computed as explained 
in the text. The next two panels show p{j) for the 3 planet resid- 
uals and the 1 planet residuals (together with a simulation thick 
solid line) . The bottom panel shows the number of HARPS sample 
pairs available for computing the autocorrelation function versus 
lag. 



is nicely accounted for by additional signals (mostly planets) 
not included in the model. 



4.6 Seven planet model 

In spite of the absence of any obvious signal from the 6 
planet residual periodogram, we computed three 7 planet 
Kepler periodogram to see if there were any new surprises. 




3.5 7.2 



28.1 38.8 53.2 91.3 



Periods 



Figure 18. A plot of eccentricity versus period for the 7 planet 
Kepler periodogram of the full HARPS data. 



Each trial used a starting period set of 7.2, 28.1, 30.8, 35.2, 
53.2, 91.3 d and a different choice for the seventh period. 
Encouragingly, all trials recovered the 6 periods found in 
the 6 planet fit which indicates they are stable features. 
Fig. [TS] shows the plot of eccentricity versus period for the 
trial which achieved the highest peak Logio [Prior x Likeli- 
hood]. The additional 3.5 d period is a weak high eccentric- 
ity feature typical of noise ( Gregory_&L_Fischcr 2010). The 
median value of the extra noise parameter s = 0.23 m s~^. 



5 MODEL SELECTION 

One of the great strengths of Bayesian analysis is the built- 
in Occam's razor. More complicated models contain larger 
numbers of parameters and thus incur a larger Occam 
penalty, which is automatically incorporated in a Bayesian 
model se lection analysis i n a quantitative fashion (see for 
example, iGregorvl (|2005al ). p. 45). The analysis yields the 
relative probability of each of the models explored. 

To compare the posterior probability of the i**" planet 
model to the 4 planet model we need to evaluate the odds 
ratio. On = p{Mi\D, I)/p{M4\D, I), the ratio of the poste- 
rior probability of model Mi to model M4. Application of 
Bayes theorem leads to. 



0»4 = 



p(M,|J) p{D\M„I) ^ p{M,\I) 
p(Mi\I) p{D\M4,I) ~ p{Ma\I) 



(8) 



where the first factor is the prior odds ratio, and the second 
factor is called the Bayes factor, Bn. The Bayes factor is 
the ratio of the marginal (global) likelihoods of the models. 
The marginal likelihood for model Mi is given by 



p{D\M„I) 



dXp{X\M,,I) xp{D\X,M„I). 



(9) 



Thus Bayesian model selection relies on the ratio of marginal 
likelihoods, not maximum likelihoods. The marginal likeli- 
hood is the weighted average of the conditional likelihood, 
weighted by the prior probability distribution of the model 
parameters. 

The marginal likelihood can be expressed as the prod- 
uct of the maxim um likelihood and the Occam penalty (see 
iGregorvl ((20053), pa-ge 48). The Bayes factor will favor the 
more complicated model only if the maximum likelihood ra- 
tio is large enough to overcome this penalty. In the simple 
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case of a single parameter with a uniform prior of width 
AX, and a centrally peaked likelihood function with char- 
acteristic width 5X, the Occam factor is « SX/ AX. If the 
data is useful then generally SX <^ AX. For a model with 
m parameters, each parameter will contribute a term to the 
overall Occam penalty. The Occam penalty depends not onf 
on the number of parameters but also on the prior range 
of each parameter (prior to the current data set, D), as sym- 
bolized in this simplified discussion by AX. If two models 
have some parameters in common then the prior ranges for 
these parameters will cancel in the calculation of the Bayes 
factor. 

To make good use of Bayesian model selection, we need 
to fully specify priors that are independent of the current 
data D. In most instances we are not particularly interested 
in the Occam factor itself, but only in the relative probabil- 
ities of the competing models as expressed by the Bayes fac- 
tors. Because the Occam factor arises automatically in the 
marginalization procedure, its effect will be present in any 
model selection calculation. Note: no Occam factors arise in 
parameter estimation problems. Parameter estimation can 
be viewed as model selection where the competing models 
have the same complexity so the Occam penalties are iden- 
tical and cancel out. 

The MCMC algorithm produces samples which are in 
proportion to the posterior probability distribution which 
is fine for parameter estimation but one needs the pro- 
portionalit y constant f o r est imating the model marginal 
likelihood. IClvde et alj l|2007t) reviewed the state of tech- 
niques for model select ion from a statistical perspective and 
iFord fc Gregorvl (|2007l ') have evaluated the performance of 
a variety of marginal likelihood estimators in the exoplanet 
context. 

Estimating the marginal likelihood is a very big chal- 
lenge for models with large numbers of parameters, e.g., our 
seven planet model has 38 parameters. In this work we em- 
ploy the nested_ restricted Monte Carlo (NRMC) method 
introduced i n iGre gorv & Fische'3 (|2010h and described in 
more detail in l Gregorvl (|2012i ') to estimate the marginal likeli- 
hoods. Monte Carlo (MC) integration can be very inefficient 
in exploring the whole prior parameter range because it ran- 
domly samples the whole volume. The fraction of the prior 
volume of parameter space containing significant probability 
rapidly declines as the number of dimensions increase. For 
example, if the fractional volume with significant probabil- 
ity is 0.1 in one dimension then in 38 dimensions the frac- 
tion might be of order 10"''*. In restricted MC integration 
(RMC) this is much less of a problem because the volume 
of parameter space sampled is greatly restricted to a region 
delineated by the outer borders of the marginal distributions 
of the parameters for the particular model. 

In RMC, most of the random samples occur close to 
the outer borders of the restricted region because the con- 
tribution to the volume of parameter space is greatest there. 
In NRMC integration, multiple boundaries are constructed 
based on credible regions ranging from 30% to 99%, as 
needed. We are then able to compute the contribution to 
the total integral from each nested interval and sum these 



^ Tlie more surprising tlie result the stronger tlie evidence re- 
quired to overcome our skepticism. 



contributions. For example, for the interval between the 30% 
and 60% credible regions, we generate random parameter 
samples within the 60% region and reject any sample that 
falls within the 30% region. Using the remaining samples we 
can compute the contribution to the NRMC integral from 
that interval. 

The left panel of Fig. \W\ shows the contributions from 
the individual intervals for 5 repeats of the NRMC eval- 
uation for the 6 planet model. Note the large range in 
parameter volume on the abscissa. The right panel shows 
the summation of the individual contributions versus the 
volume of the credible region. The credible region listed 
as 9995% is defined as follows. Let Xjjgg and Xl99 corre- 
spond to the upper and lower boundaries of the 99% cred- 
ible region, respectively, for any of the parameters. Simi- 
larly, Xj795 and Xl95 are the upper and lower boundaries of 
the 95% credible region for the parameter. Then Xygggs = 

X(799 + (Xj799 — Xj795 ) and Xl9995 — Xl99 + {Xl99 — Xl95 ) . 

Similarly, XJ79984 = Xu99 + {Xu99 - Xusi). For each credi- 
ble region interval approximately 320,000 MC samples were 
used. 

The mean value of the prior x likelihood within the 
30% credible region is a factor of 4.8 x 10^^ larger than the 
mean in the shell between the 97 and 99% credible regions. 
However, the volume of parameter space in the shell between 
the 97 and 99% credible regions is a factor of 5.9 x lO'^^ 
larger than the volume within the 30% credible region so 
the contribution from the latter to the marginal likelihood 
is neglig i ble. F or further details on the NRMC method see 
iGregorvi (|2012l 'l. 

The NRMC method is expected to underestimate the 
marginal likelihood in high dimensions and this underes- 
timate is expected to become worse the larger the num- 
ber of model para meters, i.e. increasing number of plan- 
ets (|Gregorvll20l"^ ). When we conclude, as we do, that the 
NRMC computed odds in favor of the six planet model com- 
pared to the four planet model is 137 we mean that the true 
odds is > 137. Thus the NRMC method is conservative. One 
indication of the break down of the NRMC method is the 
increased spread in the results for repeated evaluations. 

We can readily convert the Bayes factors to a Bayesian 
False Alarm Probability (FAP) which we define in equa- 
tion [TOl For example, in the context of claiming the detec- 
tion of m planets the FAPm is the probability that there are 
actually fewer than m planets, i.e., m — 1 or less. 



FAPm ~ ^ ^ (prob.of i planets) 



(10) 



If we assume a priori that all models under considera- 
tion are equally likely, then the probability of each model is 
related to the Bayes factors by 



p{PU I D, I) = 



(11) 



J4 



where A'^ is the maximum number of planets in the hypoth- 
esis space under consideration, and of course B44 = 1. For 
the purpose of computing FAPm we set N = m. Substitut- 
ing Bayes factors from Table [5] into equation 1101 gives 

p^p^ ^ (-B04 + Bi4 + B24 + + B44 + 554) ^ Q 



J4 
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Table 2. Marginal likelihood estimates, Bayes factors relative to model 4, and false alarm probabilities. The last two columns list the 
median of extra noise parameter, s, and the RMS residual. The MAP value of s appears below in parentheses. 



Model 


Periods 
(d) 


Marginal 
Likelihood 




Bayes factor 
nominal 


False Alarm 
Probability 


s 

(ms-i) 


RMS residual 
(m s-i) 


Mo 




5.69 X 10-223 




1.1 X 10-35 




3.6 
(3.6) 


3.9 


Ml 


(7.2) 


(8.57+0;«I) X 10 


-202 


1.6 X 10-" 


6.6 X 10-22 


2.3 
(2.2) 


2.6 


M2 


(7.2, 106) 


(2.37tO;«6) X 10 


-196 


6.5 X 10-9 


3.6 X 10-® 


1.8 
(1.6) 


2.3 


M3 


(7.2,28.1, 184) 


(1.19t0;2|) X 10 


-190 


2.28 X 10-2 


2.0 X 10-*^ 


1.4 
(1.3) 


2.0 


Mi 


(7.2,28.1, 106, 190) 


(5-22x2;5g) X 10 


-188 


1.0 


2.3 X 10-3 


1.0 
(0.71) 


1.8 


Ms 


(7.2,28.1,53,86.6,91) 


(3.33^J:^3) X 10 


-188 


0.64 


0.61 


0.79 
(0.70) 


1.6 


Me 


(7.2,28.1,30.8,38.8,53,91) 




186 


137 


0.012 


0.54 
(0.08) 


1.5 


Mr 


(3.5, 7.2, 28.1, 30.8, 38.8, 53, 91) 


(^■OKo.le) X 10 


-188 


0.96 


0.993 


0.23 
(0.03) 


1.4 
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Figure 19. Left panel shows the contribution of the individual nested intervals to the NRMC marginal likelihood for the 6 planet model. 
The right panel shows the integral of these contributions versus the parameter volume of the credible region. 



(12) 

For the 6 planet model we obtain a low FAP ~ IQ-^. 

Table[2]gives the NMRC Marginal likelihood estimates, 
Bayes factors and false alarm probabilities for 0, 1, 2, 3, 4, 
5, 6, and 7 planet models which are designated Mo, • • • , Mr. 
The last two columns list the median estimate of the extra 
noise parameter, s (MAP value in parentheses below), and 
the RMS residual. For each model the NRMC calculation 
was repeated 5 times and the quoted errors give the spread 
in the results, not the standard deviation. The Bayes factors 
that appear in the third column are all calculated relative 
to model 4. 

A summary of the 6 planet model parameters and their 
uncertainties are given in Table [S] The quoted value is the 
median of the marginal probability distribution for the pa- 



rameter in question and the error bars identify the bound- 
aries of the 68.3% credible region 0. The value immediately 
below in parenthesis is the MAP estimate, the value at the 
maximum of the joint posterior probability distribution. For 
the eccentricity parameter we also include the mode within 
double parentheses. It is not uncommon for the MAP esti- 
mate to fall close to the borders of the credible region. In one 



° In practice, the probability density for any parameter X is rep- 
resented by a finite list of values pi representing the probability in 
discrete intervals 5X. A simple way to compute the 68.3% credi- 
ble region, in the case of a marginal with a single peak, is to sort 
the Pi values in descending order and then sum the values until 
they approximate 68.3%, keeping track of the upper and lower 
boundaries of this region as the summation proceeds. 
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Figure 20. The marginal distribution of the 91.26 d signal (solid) 
is compared to the marginal distribution of the 30.82 d signal 
(dashed) after multiplying its period scale by a factor of 3. The 
dashed and solid vertical lines indicate the 99% credible regions 
of the two signals. 



case, the eccentricity of the fourth planet, the MAP estimate 
falls just outside the 68.3% credible region which is one rea- 
son why we prefer to quote median or mode values as well. 
The semi-major axis and M sin i values are derived from the 
mode l parameters assuming a stel lar mass of 0.31 ± 0.019 
Mq (|Anglada-Escude et al.l l2012al ) . The quoted errors on 
the semi-major axis and Msini include the uncertainty in 
the stellar mass. 



6 DISCUSSION 

We first consider whether any of the 6 signals detected is 
the result of stellar activity. D elfosse et al.i (|2oT3 ) analyzed 
several stellar activity diagnostics and found a high peak in 
one diagnostic (the full-width at half-maximum (FWHM) 
of the crosscorrelation function) with a period of ~ 105 d, 
which they interpret as the rotation period of the star. In 
our Kepler periodgrams a 106 d period was detected starting 
at the 2 planet model which transitioned to a 53 d period 
commencing with the 5 planet model. Clearly these two pe- 
riods are harmonically related. If we assume that the 106 d 
period is the star rotation period then depending on the con- 
figuration of surface activity (eg., spots) it would not be too 
surprising to detect RV variations at the second harmonic 
of the rotation frequency. This is a possible interpretation 
of the 53 d period feature in our analysis. 

In Sec. 14.51 we established that the 30.82 d period is not 
an alias of the 28.14 d period. Fig. I20l demonstrates that the 
30.82 is not a harmonic of the 91.26 d period, both of which 
show up in the 6 plan and 7 planet Kepler periodograms. 

As another test of the stability of the 6 planet model 
results, we subtracted one of the signals at a time (for the 
91, 53, and 28 d signals) from the data and carried out a 
5 planet Kepler periodogram of the modified data starting 
from an initial period set of 5, 10, 15, 20, and 25 d. In each 
case, we recovered the other 5 periods. Fig. [21] shows the 
5 planet periodogram result of the data with the 53 d sig- 
nal subtracted. The upper panel is a plot of the Logio [Prior 
X Likelihood] versus iteration for a 5 planet Kepler peri- 
odogram. The middle panel shows the values of the five pe- 
riod parameters versus iteration number. The five starting 
periods of 5, 10, 15, 20, 25 d are shown on the left hand side 
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Figure 21. The upper panel is a plot of the Logio [Prior X Like- 
lihood] versus iteration for a 5 planet Kepler periodogram of the 
HARPS data with 53 d signal subtracted. The middle shows the 
values of the five period parameters versus iteration number. The 
five starting periods of 5, 10, 15, 20, 25 d are shown on the left 
hand side of the plot at a negative iteration number. The bottom 
panel is a plot of the eccentricity parameter versus period for post 
burn-in iterations. 



of the plot at a negative iteration number. The bottom panel 
is a plot of the eccentricity versus period for post burn-in 
iterations. The marginal distributions for the 7.2, 28.1, 30.8, 
38.8, & 91 d signals are indistinguishable from those found 
for the 6 planet Kepler periodogram of the original data. 
This example also demonstrates the ability of the FMCMC 
based Kepler periodogram to carry out a blind search in 
28 parameters, including 5 period parameters which have a 
huge prior range. In fact, detection of the six signals found 
in Sec. l4.5l would be a difficult feat for any analysis method 
that relies on a conventional periodogram of the residuals 
of an n signal fit to estimate the starting period for the 
n + 1 signal search. The ability of the FMCMC algorithm 
to explore a multi-modal environment has proven to be very 
important in this analysis. 

The prospect exists that digging deeper with a more 
powerful statistical algorithm might also uncover new 
Keplerian-like spectral artifacts of the radial velocity mea- 
surement system or of the star's surface activity. There- 
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Table 3. Six planet model parameter estimates. The value immediately below in parenthesis is the MAP 
estimate. For the eccentricity parameter, the value within double parentheses is the mode. 



Parameter 


planet 1 


planet 2 


planet 3 


planet 4 


planet 5 


planet 6 


F (A) 


7.1980_ ooog 
(7.1972) 


28.1d8_ 0023 
(28.120) 


•in fi9+-04 
30.82_ Q4 

(30.82) 


q« 09+0.09 
(38.87) 


CO r,r, + 0.08 

(53.23) 


Q1 9R+0.30 

(91.37) 


A i^m s ) 


o Qc+o.ie 

•^•='"-0.17 

(3.95) 


9 91+0.21 
^•^^-0.22 
(2.32) 


1 41+0.19 
^■^^-0.20 
(1.45) 


-1 nY+O.iB 
l-'J'-0.20 

(1.33) 


1.40_o 21 
(1.55) 


(1.76) 


e 


(0.068) 
((0.068)) 


0.08d_p gg^ 

(0.123) 
((0.028)) 


(0.22) 
((0.115)) 


n ^c:+0.21 

(0.58) 
((0.48)) 


n 4-1+0-12 
'J-^i-o.io 
(0.48) 
((0.45)) 


n qc+O.lO 

(0.34) 
((0.35)) 


^ (deg) 


342^37 
(338) 


269l«« 
(266) 


209l« 
(220) 


284^31 
(297) 


3051^3 
(296) 


219t^o 
(217) 


a (AU) 


0.0491:001 
(0.0494) 


123+°°^ 
(0.123) 


0.130100^ 
(0.130) 


U.10Z_ 003 

(0.152) 


n -IS7+0.004 

(0.187) 


n 9«o+0.006 

(0.269) 


Msini (Mq) 


■5 4+0.3 
(5.44) 


^•^-0.5 

(5.01) 


3 1+0* 
(3.17) 


2 4+0-4 
(2.62) 


3 4+0.5 
(3.67) 


5 4+0.5 
(5.34) 


Periastron 
passage 


4495.6to;? 
(4503.0) 


4472^5 
(4499) 


4468^1 
(4503) 


4480^2 
(4480) 


4447^2 
(4499) 


4387^4 
(4477) 



(JD - 2,450,000) 



fore, N-body simulations will be required to determine which 
combination of the signals are consistent with a stable plan- 
etary system. 

In the rest of this discusion we speculate on the pos- 
sibility that some or all of the remaining 5 signals de- 
tected are planetary in origin. One surprise is the close- 
ness of the 28 d and 30.8 d orbits whose semi-major axes 
differ by only 0.007 AU. Closest approach would occur ev- 
ery 1/(1/28.138 - 1/30.820) = 323 d. There is precedence 
for systems with very close separations. One of the Kepler 
mission's many interesting findings is that many stars host 
multiple planets in surprisingly close orbits. Examples in- 
clude : Kepler 36b & c (M = 4.3 & 7.8 M0, ICarter et all 
|2012| ) with a minimum separati on of 0.014 AU, Kepl er 42b 
& d (M ^ 2.7 & < 0.9 M^. iMuirhead et aLllioij ) with 
a minimum separation of .0038 AU, and KOI-5 5b & c 
(M ~ 0.44 & ~ 0.66 M^. ICharpinet et all [20]! ) with a 
minimum separation of 0.0016 AU. The masses of the par- 
ent star for Kepler 36, Kepler 42 and KOI-55 are 1.07, 0.13 
and 0.496 Mq, respectively. Kepler 36b & c both show tran- 
sit timing variations consistent with a gravitational interac- 
tion between the two planets while no such timing variations 
were detected for the three Kepler 42 planets. 



6.1 Multiple planets in the habitable zone? 

By definition, the habitable zone around a host star is the 
region where liquid water can be stable on the surface of a 
rocky planet l|Huanelll959l . iKasting et aL 19931). In the cur - 
rent absence of observational constraints, Selsis et al.l (|2007l ) 
chose to assess the habitable zone for planets with as few as- 
sumptions as possible on their physical and chemical nature. 
Only two conditions were assumed although they may derive 



from complex geophysical properties. Below we summarize 
these assumptions and some of the key concepts. 

1) The amount of superficial water must be large enough so 
that the surface can host liquid water for any temperature 
between the temperature at the triple point of water, 273 K, 
and the critical temperature of water, 647 K. This condition 
implies that the water reservoir produces a surface pressure 
higher than 220 bars when fully vaporized. With an Earth 
gravity, this corresponds to a 2.2 km layer of water, slightly 
lower than the mean depth of Earth oceans of 2.7 km. For a 
gravity twice that of Earth, this pressure corresponds to half 
this depth. Planets with less water may still be habitable, 
but their HZ may be somewhat narrower because liquid wa- 
ter would disappear at a lower surface temperature. 

2) Atmospheric CO2 must accumulate in a planet's atmo- 
sphere whenever the mean s urface temper ature falls below 
the freezing point of water. IWalker et al.l (1981) proposed 
that Earth's long-term climate was buffered by a negative 
feedback mechanism involving atmospheric CO2 levels and 
the dependence of silicate weathering rates on climate. In 
the carbonate-silicate cycle, the CO2 emitted by volcanoes 
is dissolved in rain water forming carbonic acid (H2CO3) 
which weathers silicate rocks. The products of silicate weath- 
ering, calcium (Ca"*"^), bicarbonate (HCO;^) ions and dis- 
solved silica (3102), are transported by streams and rivers 
to the ocean. There, organisms, such as foraminifera, di- 
atoms and radiolarians, use the products to make shells of 
calcium carbonate (CaC03). Most of the shells redissolve, 
but a fraction of them survive and are buried in sediments 
on the seafloor. The combination of silicate weathering plus 
carbonate precipitation can be represented chemically by 
CO2 -I- CaSi03 -> CaCOs -I- Si02. The seabed is eventu- 
ally subducted where the heat pushes the carbonate-silicate 
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cycle reversible reaction in the opposite direction, eventually 
releasing CO2 through volcanism. 

This cycle replenishes all the CO2 in the combined 
atmosphere-ocean system on a timescale of approximately 
half a million years. It is thus too slow to counteract human- 
induced global warming, but fast enough to have a dominat- 
ing effect on the billion-yea r timescale of planetary evolution 
l|Kasting fc Catling \\20m . Weathering rates increase both 
because of the direct effect of temperature on chemical re- 
action rates and because evaporation (and, hence, precipi- 
tation) rates increase as the surface temperature increases. 
If we increase the orbital radius, the surface temperature 
falls and the weathering decreases allowing atmospheric CO2 
to build up leading to a stabilization of the temperature 
through the greenhouse effect. This negative feedback in the 
carbonate-silicate cycle stabilizes the long-term surface tem- 
perature and the amount of CO2 in the atmosphere of the 
Earth. Such an assumption implies that the planet is ge- 
ologically active and continuously outgassing CO2. It also 
implies that carbonates form in the presence of surface liq- 
uid water, which may require continental weathering. 

At a sufficiently large orbital radius the planet gets cold 
enough that water vapor disappears followed by the freezing 
of CO2 which results in the permanent loss of both green- 
house gases (positive feedback) . With no atm ospheric CO2, 
or with a fixed CO2 level as in iHartl l{197^ . the HZ could 
be ~ 10 times narrower. In the absence of a greenhouse gas 
like CO2 and H2O, the present Earth would be frozen. 

With decreasing orbital radius, a positive feedback 
eventually takes over in which too much water accumulate 
in the atmosphere causing the temperature to increase. At a 
sufficiently high temperature water stops condensing. With 
no rain the weathering ceases and CO2 builds up, further 
increasing the temperature. Eventually the oceans evapo- 
rate completely. Water in the high atmosphere is dissociated 
into hydrogen and oxygen, the hydrogen escapes to space, 
the oxygen combines with minerals and all the water disap- 
pears. 

According to ISelsis et al.l (|2007h , planets with masses 
outside the 0.5 - 10 M0 range cannot host liquid surface 
water. Planets under the lower end of this range have too 
weak a gravity to retain a sufficiently dense atmosphere, and 
those above the upper end accrete a massive He-H envelope. 
In either case, the pressure at the surface is incompatible 
with liquid water. The 10 upper limit is somewhat fuzzy, 
since planets in the 3-10 range can have very different 
densities. 

Based on t he above two major assumptions, 
l|Selsis et all |2007| ) derives boundaries for the HZ as 
show in Fig. [21] for a range of stellar masses. Assuming 
a planetary origin for the Keplerian-like signals reported 
here, the 28.1, 30.8, and 38.8 d signals (labeled c, d, and e) 
translate to orbits in the centre of the HZ. The semi-major 
axis of a 91.3 d planet (f) would lie just within the outer- 
most extreme edge of the HZ, although its eccentric orbit 
will result in it spending the majority of its time outside 
the HZ. The white object in the diagram corresponds to the 
orbital parameters implied by a 53 d Keplerian signal. As 
explained above, it is likely that the star's surface activity 
is responsible for this signal as it is the second harmonic of 
what is believed to be the star's rotation period. 




0.1 1.0 10.0 

Distance from star (AU) 

Figure 22. The darker area is the orbital region that remains 
continuously habitable during at least 5 Gyr as a function of the 
stellar mass l lSelsis et al.ll200^ . The light grey region gives the 
theoretical inner (runaway greenhouse) and outer limits with 50% 
cloudiness, with H2O and CO2 clouds, respectively. The dotted 
boundaries correspond to the extreme theoretical limits, found 
with a 100% cloud cover. The dashed line indicates the distance 
at which a 1 Mgj planet on a circular orbit becomes tidally locked 
in less than 1 Gyr. 



7 CONCLUSIONS 

A Bayesian r e-analysis of published HARPS precision radial 
velocity data iDelfosse et al.l l|2012l ) for Gl 667C was carried 
out with a multi-planet Kepler periodogram (from to 7 
planets) based on our fusion Markov chain Monte Carlo al- 
gorithm. In all cases the analysis included an unknown pa- 
rameterized stellar jitter noise term and an unknown slope 
parameter to deal with the linear drift resulting from the 
orbital motion of Gl 667C relative to the center of mass of 
the Gl 667ABC system. The most probable number of sig- 
nals detected is 6 with a Bayesian false alarm probability of 
0.012. The residuals are shown to be consistent with white 
noise. The Keplerian signals detected include the 7.2 and 
28.1 d (planets b and c) periods reported previously plus 
periods of 30.8 (d), 38.8 (e), 53.2, and 91.3 (f) d. The 53.2 
d period appears to correspond to the second harmonic of 
the stellar rotation period and is likely the result of surface 
activity (eg., spots). The same set of periods were also de- 
tected in a subset of the data consisting of the first 143 data 
points. 

The existence of the additonal Keplerian-like signals 
suggest the possibilty of further planets, two of which (d 
and e) could join Gl 667Cc in the central region of the hab- 
itable zone. Msini values corresponding to signals b, c, d, 
e, and f are ~ 5.4, 4.8, 3.1, 2.4, and 5.4 M0, respectively. 
It is also possible that digging deeper with a more powerful 
statistical algorithm might have uncovered new spectral ar- 
tifacts of the radial velocity measurement system or of the 
star's surface activity. N-body simulations are required to 
determine which of these signals are consistent with a sta- 
ble planetary system. 
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